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Abstract. This paper constructs dynamical models and estimation algorithms for the con- 
centration of target molecules in a fluid flow using an array of novel biosensors. Each biosensor 
is constructed out of protein molecules embedded in a synthetic cell membrane. The concentra- 
tion evolves according to an advection-diffusion partial differential equation which is coupled with 
chemical reaction equations on the biosensor surface. By using averaging theory methods and the 
divergence theorem, an approximate model is constructed that describes the asymptotic behaviour 
of the concentration as a system of ordinary differential equations. The estimate of target molecules 
is then obtained by solving a nonlinear least squares problem. It is shown that the estimator is 
strongly consistent and asymptotically normal. An explicit expression is obtained for the asymptotic 
variance of the estimation error. As an example, the results are illustrated for a novel biosensor built 
out of protein molecules. 
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1. Introduction. Estimating the concentration of target molecules in a fluid 
flow over multiple biosensors is a challenging problem due to two non-standard fea- 
tures. Firstly, it is a parameter estimation problem of an advection-diffusion partial 
differential equation (PDE) which is coupled with Dirichlct and Neumann boundary 
conditions and cannot be solved analytically. Secondly, the measurement process af- 
fects the system state since each biosensor grabs target molecules and changes the 
concentration in the flow. This is unusual since in most statistical inference problems, 
observation does not change the system state. The main results of this paper are 
briefly stated as follows: 

1. An advection-diffusion PDE model is constructed to model the variations 
of the concentration of target molecules that flow past a linear array of biosensors. 



To facilitate estimation of the concentration of target molecules, Theorem 3.1 devel 



ops an approximation method to describe the dynamics of the problem by a system 



of ordinary differential equations (ODEs). In Theorem 3.1 the multi-compartment 
ODE model is derived by exploiting the multiple time-scale behaviour of the system, 
together with perturbation methods and the divergence theorem. 

2. A novel biosensor constructed out of protein molecules is used as an actual 
example to illustrate our results. The development of this biosensor was first published 
in Nature \1 . The biosensor incorporates ion channels into a tethered lipid bilayer 
membrane where the conductance of the channels is switched by the recognition event. 
The PDE model for this biosensor is specified and solved numerically using the Comsol 
multi-physics finite element analysis software. We show how the PDE model and ODE 
approximations can satisfactorily model this novel biosensor. 

3. The estimation of target molecule concentration is posed as a parameter es- 
timation problem in terms of the derived ODE model. The estimate is computed nu- 
merically for the novel biosensor via nonlinear least squares method. The asymptotic 
behaviour of the estimator is analyzed. It is shown that the estimator is asymptoti- 
cally unbiased and normal and its asymptotic variance is derived. According to the 
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expression for the variance, the achievable improvement in the estimate based on the 
number of biosensors is evaluated using the results from the ODE model. 
Inverse problems in fluid mechanics deal with estimating an unknown coefficient or 
function in the initial or boundary condition for a parabolic PDE (see [2-6 ). In 
our case, the problem is estimating the boundary condition of a parabolic PDE. The 
model-based state or parameter estimation of distributed systems (infinite dimensional 
systems) based on a distributed-parameter description is quite complex. To address 
this problem, the system description is converted from a distributed-parameter into 
a lumped-parameter form. This conversion can be achieved by methods for solving 
partial differential equations, such as finite-difference method [7], the finite-element 
method, modal analysis [8] and spectral method |9j. 

The estimation problem in this paper can also be viewed as a source determination 
in a fluid system. In [10] , the strength (emission rate) of a contaminant source is 
estimated using a fixed network of concentration measurements and a Lagrangian 



trajectory model. In 11 , a backward-time Lagrangian stochastic model is used to 
estimate the emission rate of a surface area source. In [12] , Baysian inference is 
applied to solve a chemical source determination problem where the posterior joint 
distribution of location, intensity, and temporal properties of a point source is obtained 
by a Markov chain Monte Carlo (MCMC) method. The authors develop a dual 
problem for the advection-diffusion equation using adjoint dispersion equations which 
requires significantly less amount of computations comparing to resolving the main 
equation for every source term. 

In the above works, the PDE has standard non-coupled initial and boundary con- 
ditions where measurements do not affect the system state whereas our problem has 
unconventional boundary conditions. In order to convert the system description from 
a distributed-parameter into a lumped-parameter form, an approximation method is 
used which is based on the two-compartment model |13| . The two-compartment model 
is used in modeling a variety of binding experiments influenced by diffusion and mass 



transport. For example, in 14 and 15 , this model is used to study and characterize 
the kinetic properties of biomolecular interactions in optical biosensors. In this work, 
the two-compartment model is extended and developed to a new multi-compartment 
model to describe the reaction-diffusion experiment in a flow chamber over a linear 
array of multiple biosensors. With this method, the PDE model is approximated by 
a system of ordinary differential equations. 

In Secj2] the PDE model for a general reactive surface is described. It is followed 
by the derivation of the multi-compartment ODE model in Secj3] The asymptotic 
properties of the least squares estimator for the concentration of target molecules 
is studies in Sec]4] using the multi-compartment model of Sec[3] Sec[5] presents the 
results for a protein-based biosensor. 

2. PDE model for the fluid flow. The aim is to estimate the concentration of 
target molecules in a fluid system where the dynamics are described by an advection- 
diffusion PDE model. 

Consider a flow chamber with a rectangular cross section where a flow of tar- 
get molecules flows past multiple surface-based biosensors along the length of the 
chamber. There are N identical biosensors which form a linear array along the flow 
direction on the surface of the chamber floor. We introduce three-dimensional Carte- 
sian coordinates (x, y, z) with the y-axis along the flow direction and the z-axis along 
the height of the flow chamber and perpendicular to the surface of the biosensors. 
Biosensor i, for i = 1, 2, . . . , N, is located in the range Vi,2\ along the y-axis and 
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Fig. 2.1. An equally spaced linear array of three biosensors in a rectangular flow chamber. 
The fluid containing target molecules enters from the left side. Th e con centration of target 
molecules at the inlet is A* as expressed in the boundary condition (2.5\). 



[0, w] along the s-axis. The inlet of the flow chamber lies in the x — z plane. The 
system is symmetric about the x-axis since the ratio of the height to the width of 
the flow chamber is selected to be less than 1/20 [16]. The dimensions of the flow 
chamber and biosensors are: 

Flow chamber: Height = h, Length = I, Width = uj, (2.1) 
Biosensors: Length = L, Spacing = d, 

Biosensor i is located in the range y € [2/1,1,2/1,2] , < 2/1,1,2/1,2 < I- 



A flow chamber with N = 3 biosensors is illustrated in Fig |2.1| When target molecules 
in the solution arrive at the biosensors, chemical reactions are initiated which result 
in a change in impedance that is translated to change in the measured current. 

Below, an advection-diffusion PDE is used to describe the spatio-temporal evo- 
lution of concentration of target molecules in the flow chamber. It is coupled with a 
set of ODEs on the boundary which describes the adsorption of target molecules on 
the biosensors as a result of chemical reactions. 

Fluid flow dynamics: The concentration of target molecules in the flow cham- 



ber (2.1), denoted by A(t,y,z) is governed by an advection-diffusion PDE 17 



OA fd 2 A d 2 A\ . BA , n n , n 

M = \W + B*)- V MW Vei °' l) ' ZeM - (2 ' 2) 

Here 7 is the diffusion constant of the target molecule and v(z) is the flow velocity in 
y direction. The flow is assumed to be laminar and fully developed with a parabolic 
velocity profile defined by 

v(z) =Av{z/h){l~ z/h), (2.3) 

where v is the maximum velocity [15] . Initially, the flow chamber is empty and the 
concentration inside the flow chamber is zero. So the initial condition is written as 

A(t = 0,y,z) = 0, 2/6(0,0, z 6(0, ft). (2.4) 

To complete the modeling of the fluid flow, it is necessary to specify the following 
Dirichlet and Neumann boundary conditions: The concentration at the inlet of the 
flow chamber is constant during the estimation process and equal to A* . At the outlet 
of the flow chamber, there is no diffusive flux. There is insulation at the ceiling of 
the flow chamber at z = h. The gaps between the sensors on the floor of the flow 
chamber, at z = 0, are also insulated. These boundary conditions are described as 



8 A 

A(t iy = 0,z) = A*, — 
oy 



= d A 

y=i dz 



= d A 

z=h dz 



= 0. (2.5) 



On the surface of each biosensor, the adsorption flux of target molecules is equal to 
the rate of consuming target molecules by the reactions. 
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Chemical dynamics: Assuming that biosensor i, i € {1,2, . . . , N}, is located 
in y £ [yi,i, Ui,2\ on y-axis, the corresponding boundary condition is expressed as 



17 



dA 

az 



= R(A(t,y,z = 0),u l (t,y)), 



(2.6) 



z=0,v€[Vi,i,Vi,2] 



where the vector Uj(t, y) contains the values of concentration of chemical species at 
time t and location y on biosensor i. R(A, Uj) is the rate of adsorption of target 
molecules per unit area on the biosensor surface. The rate of adsorption at each point 
is a function of the concentration of target molecules and chemical species on that 
point on the biosensor. The dynamics of the chemical species of biosensor i at location 
y are described by a system of ODEs as 



duj(t,y) 
dt 



G(ui(t,y),A(t,y,0)), t > **, Ui(U,y) = uo, for y G [^,1,^,2] • (2.7) 



Here, G(-) is a function which is described by the rate law of reactions on the biosensor, 
ij is the time instant at which the flow reaches biosensor i and the biosensor starts 
responding. The rate of change of Uj(t, y) depends on the concentration A{t, y,z = 0) 
of target molecules on the biosensor. The constant uq is the initial concentration of 
chemicals on each biosensor. 

Aim: The aim is to estimate the concentration A* at the inlet of the flow cham- 
ber which appears in the boundary condition (2.5|. After describing the biosensor 



array model, statistical estimation algorithms are given in Sec(3]and Sec|4]to estimate 
A* given noisy measurements from the biosensors. 

Measurement equation: Finally, the measurement equation is specified. The 
response of biosensor i at time t is denoted by gi{A\,t) where A\ refers to the con- 
centration at the inlet. The response gi(Ai,i) can be described as a function of Qj(t) 
which denotes the surface average of the concentration vector Uj(t, y) on biosensor i; 



gi{Ax,t) = F(ui(t)), Ui(t) 



1 

Ui,2 - y%,\ 



Vi,2 



Ui(y,t) dy. 



(2.8) 



ViA 



Here, F(-) is the transducer function which translates the concentration quantities 
on the biosensor to a corresponding electrical signal. In Sec(5j we give a specific 
example of an actual biosensor where F(-) models the conductance of the biosensor. 
Considering the PDE model (2.2|-(2.7|, the measurement taken at biosensor i at time 
t l ' k , denoted by mf, is 



= g i (A*,t i < k ) + nl ie{l,2,...,N}, k e {1, 2, . . . , S}. 



(2.9) 



Recall A* is the value of the concentration Ai at the inlet and is the corresponding 



measurement noise. In (2.9), S is the number of measurement samples taken at each 
biosensor. The noise samples vq for i = 1, . . . , N and k = 1, . . . , S are independent 
normally distributed with zero mean and finite variance a 2 . In (2.9 1, gt(A* , t l > ) is an 
implicit function of the concentration A* through the PDE model (2.2)-(|2.7|). 



3. Multi-compartment model approximation. Given the measurement equa- 
tion of (2.9) and the PDE model of Secj2]defined by (2.2)-(2.7), the aim is to estimate 
the concentration A* at the boundary i n (|2.5[ ). The PDE is coupled with a set of 
ODEs through the boundary condition ( |2.6[ ) and cannot be solved analytically. To 
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Fig. 2.2. Four- compartment model for two biosensors. 



estimate A* in (2.5), in this section, a multi-compartment ODE model is introduced 
that approximates the PDE by a system of ODEs. 

The multi-compartment model is an extension of the existing two-compartment 
model to a new model for mass transport-binding experiments on a linear array of 
multiple biosensors. Its derivation is based on the multiple time-scale behaviour of 
the system [18] and the divergence theorem. First, the two-compartment model is 
reviewed in Sec|3.1| Then, the multi-compartment model is then derived in Sec|3.2| 



3.1. Review of the two-compartment model. The two-compartment model 
consists of a set of coupled ordinary differential equations which are used to analyze a 
variety of binding experiments influenced by mass transport 



15 . It is derived based 



on the multiple time-scale behaviour of the system and models the slow response 
of the system. When the intrinsic reaction rates are comparable to or faster than 
the rate of transport of molecules to the reactive surface, a depleted region forms 
on top of the reactive surface where the amount of target molecules is slowly varying 
comparing to the concentration in the bulk. Because of this two time-scale behaviour, 
the flow chamber on top of the biosensor is divided vertically into two compartments, 
as shown in Fig |2.2[ in order to consider their dynamics separately. This model ignores 
the brief transitions before the bulk (outer compartment) concentration falls or rises 
to the concentration A* at the inlet. The concentration of target molecules within 
each compartment is assumed to be spatially homogeneous. The concentration in the 
outer compartment, denoted by A±, is equal to the concentration A* at the inlet of 
the flow chamber. The dynamics of the average concentration of target molecules in 
the inner compartment, denoted by <xi(t), is described by 



13 



ho 



where 



dt 



-i?(ai(i),u 1 (i)) + 7 I (A 1 -a 1 (i)), ,i>0, oi(0)=0, (3.1) 

tlQ 



h = 



1 

1.464 



jhL' 
v 



1/3 



(3.2) 



is the height of the inner compartment. In (3.2 1, 7 is the diffusion constant, v is the 
maximum velocity in the velocity profile (2.3), and R{-,-) is defined in (2.6). The 



average concentration of chemicals on the biosensor, denoted by Ui(t), is governed by 
dui(i) 



dt 



G(iii(t),ai(t)), t>0, ui(0) 



u . 



(3.3) 



3.2. The multi-compartment model. In this section, the PDE model (2.2) 



( 2.7 ) is approximated by a system of ordinary differential equations. The flow chamber 
is partitioned into a series of two-compartment blocks above biosensors which are con- 
nected by middle compartments as shown in Fig |2.2| As shown in Theorem [3T] below, 
the response of each biosensor can be described by an individual two-compartment 
model with a different concentration in its outer compartment. We apply the di- 



vergence theorem to the advection-diffusion PDE of (2.2) in the outer compartment 
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associated with each biosensor to find the concentration at the inlet of the next biosen- 
sor. By exploiting the different time-scale dynamics of the concentration in the flow 
chamber and applying some perturbation methods, a two-compartment model for the 
next biosensor can be derived. The concentration in the outer compartments of con- 



secutive biosensors are related by (3.7) in Theorem 3.1 For the derivation of the 

multi-compartment model, the following assumptions are required: 

(1) The concentration A(t, y, z) in the flow chamber is increasing and concave in t for 



y < tv where v is the maximum velocity in the velocity profile (2.3). 

(2) The response of biosensor i commences at time t% where 

U = ^, i = l,2,...,N. (3.4) 

v 

Thus, ti is the time instant at which the flow reaches the far end of biosensor i at 
V — Vi,2 in the flow chamber. 

(3) The spacing of biosensor i is sufficiently small such that 7 1 ^ 3 (^ — ti-i) = 0(7 1 / 3 ) 
as 7 — > for i € {1, 2, ... , N}. Here, 0(.) denotes the Landau big-O. 

Assumption (1) models the physical reality that the concentration A(t, y, z) at each 
point is bounded and eventually asymptotes at the concentration A*. Therefore, 
A(t,y, z) is concave in time after a certain time instant. Assumption (2) models the 
two-time scale behaviour of the flow chamber: the concentration in the outer compart- 
ment of each biosensor evolves rapidly compared to that in the inner compartment 
(vicinity of the surface of the biosensor). Assumption (3) reflects the fact that the 
spacing of the biosensors should be sufficiently small such that each biosensor is af- 
fected by the depletion region that is generated by the previous one. Otherwise, the 
biosensors have identical responses. The above assumptions are justified by numer- 



ous experimental studies of the biosensor, see 17 . The following multi-compartment 
characterization is the main result of this section. 

Theorem 3.1. Consider a flow of target molecules over an equally spaced linear 



array of N identical biosensors in the flow chamber (2.1). Suppose the concentration 
of target molecules at the inlet of the flow chamber is a constant denoted by A* . The 
concentration of target molecules Ait, y, z) and chemical species are described by the 



PDE model (2.2)-(2.1). Under assumptions (1) and (2), as 7 0, there exists a 
time instant t* such that for t G (tj,t*), the dynamics of the average of the surface 
concentration of chemical species on biosensor i, denoted by Ui(t), satisfies 

hQ dagl = _7 (A . _ W) _ R( - i{th ^ {t)) + 0(7 4/ 3)j t e {UX) 
cit riQ 

dv ^ = G(u t (t), at® + o( 7 2 )), * e A 

a.i{U) — 0, Ui(ti)=u , i=l,...,N. (3.5) 



Recall"/ in (3.5) denotes the diffusion constant. The concentration in the flow chamber 



for z £ (ho, h — ho) can be expressed as 

A(t,D,z) = Ai + 0{i), y € (t/i-i,2,2/i,2), z e (h ,h- h ), t € (t,*,t*), (3.6) 
where ho is defined in $3.2\). Ai is a constant obtained by the following recursion: 



A i = aA i - 1 , i = 2,...,N, a= [ 1 - ) , A l= A*. (3.7) 

zhovn ' 
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Proof. The proof is in Appendix [A] □ 

The implication of the above theorem is that the multi-compartment model pro- 
vides an accurate description of the dynamics of the concentration of target molecules. 



In (3.5 1, Theorem 3.1 ho is the height of the inner compartment above each biosensor 
and CLiit) denotes the spatial average of concentration in the inner compartment of 
biosensor i. The non-negative constant denotes the concentration in the outer 



compartment of biosensor i. The ODEs of the two-compartment model in (3.1) and 



(3.3) are generalized for biosensor i in (3.5) in this model 



4. Asymptotic analysis of the least squares estimator of the concen- 
tration. With the above multi-compartment characterization of the concentration 
of target molecules, this section deals with estimating the initial concentration A*. 
The estimation is formulated as a least squares problem for the multi-compartment 
model. Then the asymptotic behaviour of the estimator in terms of consistency and 
asymptotic normality is investigated. Also an approximate formula for the variance 
of the finite-sample estimator is obtained. This allows us to evaluate qualitatively 
how the variance of the estimate varies with the number of biosensors. 

The concentration A\ — A* at the inlet of the flow chamber is estimated using 
non-linear regression. The estimate denoted by A\ is defined as 

N S 

ij = arg min S^VVK - ^(Ai^ k )f (4.1) 
where N and S refer to the number of biosensors and the number of time samples, 



respectively. In (4.1), m k refers to rrii(t hk ) which, according t o (|2.9 l), is the measure- 
ment of biosensor i taken at time t l,k . <fr(Ai, t*'*) defined in (|2.8[), is the respo nse o f 



biosensor i at time t l ' k when the concentration at the inlet is A\. Regarding (2 
gi(A\, t l ' k ) is a function of \ii(t). Thus, the estimate is the solution of the optimization 
problem ( |4.1[ ) together with ( |3.5[ ). 

The following definition, assumptions and theorems establish strong consistency 
and asymptotic normality of the least squares estimator for the concentration A* . 

Definition 4.1. Let f — (f t ) and g = (g t ) be two sequences of real valued 
functions on and h be a function on x 0. If ft -1 X^tLi ft{ a )9t{ft), as n — > oo, 
converges uniformly to h(a,f3) for all a and j3 in 0, h— [f,g] is called the tail cross 
product of f and g. 

Assume that: 

(a) a sequence of real valued N x 1 vectors y t has the structure y t = ft(#o) + e t 
for t = 1,2,3,... where the elements of the vector function f t , denoted by f^t for 
i = 1, 2, . . . , N, are known continuous functions on a compact subset of a Euclidean 
space and the vectors e t for t = 1,2,3,... are independent identically distributed 
with zero mean and finite covariance matrix <7 2 //v > 0. Here, In denotes the N x N 
identity matrix. 



(b) considering Definition 4.1 the tail cross product of /, = (/^t), denoted by 
[fi, fi], for i = 1,2,...,N exists. Besides, Q{6) = lim^ n" 1 \ft(9 ) - f t (0)| 2 
has a unique minimum at 9 = 9q. Here, |.| 2 denotes the L-2 vector norm. 

Any vector n in which minimizes Q n (9) — nT 1 X^t=i |y* — ^t{9)\ , is a least 
squares estimate of 9q based on the first n values of y t . The following theorem 
establishes strong consistency of 9 n . 

Theorem 4.2. Suppose that (9 n ) is a sequence of least squares estimators. Under 
assumptions (a) and (b), n is a strongly consistent estimator of9g. 
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Proof. The proof for a scalar valued function / is given in 19 . Extending the 
proof to a scalar valued function is straightforward. □ 

To establish the asymptotic normality of a sequence of least squares estimators 
of a scalar parameter, we need the derivatives f' it {9) — (d / d0)fi.t(9) and f" t {0) = 
(d 2 /d9 2 )f ht (9) for i = 1,2,..., AT and t = 1,2,3,.'... Suppose that: 

(c) the derivatives f[ t and f" t for i = 1, 2, . . . , N exist and are continuous on 
and that all tail cross products [fi, //], [/*,/"], and [/-, /"] for i = 1, 2, . . . , N exist. 

(d) for each 9 in 9, a(8) is defined as a(6) = lim, woo n" 1 J2t=i J2iLi (/i,t(^)) 2 - 
The true parameter 9q is an interior point of and a(#o) is not zero. 

The following theorem provides conditions for the asymptotic normality of a se- 
quence of least squares estimators. 

Theorem 4.3. Suppose that (9 n ) is a sequence of least squares estimators of a 
scalar parameter 9q. Under assumptions (a) through (d), 9 n — 9q is asymptotically 
normal, i.e. ^fn (§ n - O ) -> N(0, a 2 a(9 )~ 1 ). 

Proof. The proof for a scalar function / is given in |19| . Extending the proof to 
vector valued functions is straightforward. □ 

Using Theorem |4.2| and Theore m |4 .3j strong consistency and a symp totic normal- 



below. 



ity of the estimator A\, defined in (4.1 1, is established in Theorem 

Theorem 4.4. Consider the observation model (2.9), dynamics [3.5), the re- 



4.4 



lation of the biosensor response with the concentration of species in (2.S) and the 



recursion (3.1). Assume that the noise samples in (2.9) are independent identically 
distributed with zero mean and finite variance a 2 . Then the estimate A\, defined in 



(4-D, has the following asymptotic properties: 

1. A\ is strongly consistent as the time sample size S in (4-D grows. 

2. The estimation error A\ — A* is asymptotically normal as S — > oo; 



N(0, 



1 



N 



r= lim TV 



, 2i-2 



fc=l i=l 



dg(a l ~ 1 A*,t l ' k - U) 
dA 



(4.2) 



where g(A,t), defined in (2 



is the response of each biosensor when the concentra- 
tion in its outer compartment is A. dg(a l ~ 1 A* , t l ' k — ti)jdA is the value ofdg(A, t)/dA 
at A = a 1 ' 1 A* and t = V" k - U. 



a A ant 
Proof. The proof is given in Appendix [5] □ 
Corollary 4.5. Consider the observation model 



2JA), dynamics 



relation of the biosensor response with the concentration of species in (2.8) and the 



recursion (3.7). Assume that the noise samples in (2.9) are independent identically 
distributed with zero mean and finite variance a 2 



3.5), the 



Suppose that a 



i 

S,N 



variance of the finite sample estimator A±, defined in (4-D for finite S 



N 



2 

a S,N 



where §^{ct l 1 A* 



k=l 



dg(a l ~ 1 A* ,t l ' k -t t ) 
dA 



n 2 



denotes the 
Then, 

(4.3) 



t> 



-ti) is the value ofdg(A,t)/dA at A — a 1 X A* andt = t l ' k —ti 



The proof of Corollary |4.5| follows straightforwardly from Theorem 4.4 Corol 



lary |4.5| shows how the variance of the estimate of the concentration A* varies with 
the number of biosensors. In the experiments involving the protein-based biosensor 
described in Sec[5j the approximation (4.3) is used to explain how the estimation 
variance varies with the number of biosensors. 
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5. Case-study: Ion channel biosensor. In this section, the multi-compartment 
model of Sec |3.2| is evaluated for a protein-based biosensor, namely the ion channel 
switched (ICS) biosensor that was constructed and described in (20l. ICS biosensor 
is a generic biosensor that can detect low molecular weight drugs, large proteins and 
micro-organisms [20] with low concentrations as low as 10 fMolar [l7j. This biosen- 
sor incorporates artificial ion channels in a lipid bilayer. The PDE model of Sec(2] 
is specified for this biosensor by describing the corresponding chemical reactions and 



measurement equation in Sec 5.1 In Sec 5.2 the multi-compartment model (13. 5| is 



applied to an ICS biosensor array. Comparison between its response with the response 
of the PDE shows that the multi-compartment model describes the system accurately. 
Finally, in Sec|5.3[ numerical examples are given for estimating the concentration A* . 



5.1. Dynamics of the flow on (Ion Channel Switch) ICS biosensor. In 

this section, the operation of the ICS biosensor is outlined. Then, the PDE model of 
Secj2]is constructed for a linear array of ICS biosensors. 

Chemical Dynamics: Recall Secj3]gave a generic description of chemical reac- 



tions in (2.7). Here, the specific chemical dynamics on ICS biosensor are described. 



More details on the construction and operation of this biosensor can be found in 17 
and [20] . To specify chemical dynamics, we first outline briefly the structure and oper- 
ation of the ICS biosensor. The ICS biosensor is a surface based biosensor comprised 
of a lipid bilayer where ion channels are infused. The inner lipid layer is tethered to 
a gold substrate. The ion channels within this layer are tethered whereas the ones in 
the outer layer diffuse freely. The flow of ions through a channel only occurs when a 
mobile channel in the outer layer aligns to a fixed channel in the inner layer to form 
a conducting dimer. The arrival of target molecule cross-links antibodies attached to 
the mobile outer layer channels, to those attached to tethered lipids. This anchors 
the channels distant, on average, from their inner layer partners. The expected num- 
ber of dimers is thus decreased. The conductance of the biosensor is proportional 
to the concentration of the dimers. Therefore the arrival of target molecules results 
in decreasing the biosensor admittance. There are eight chemical species on the ICS 
biosensor (17) . Therefore, the vector of concentration of chemical species u has eight 
elements. The primary species include binding site b with concentration B, free mov- 
ing ion channel c with concentration C, tethered ion channel s with concentration 
S, and dimer d with concentration D. Initially, free moving ion channels c, tethered 
channel s, dimers d are in equilibrium through a reversible chemical reaction. The ar- 
rival of target molecules initiates six other reactions and the equilibrium shifts towards 
decreasing the dimer concentration. The target molecule binds to the primary species 
to form complexes w, x, y, and z with concentrations W, X, Y , and Z according to 
the following chemical reactions 

a + b^l^w a + c^l 2 o x w + c^hy x+b^^y , . 

c + s d a + d ^ z x + s z I • J 

Here, fj and rj for j — 1, 2, . . . , 7, respectively denote the forward and backward 



reaction rate constants. The corresponding rate equations for the reactions (5.1) are 



R 1 = f 1 AB-r 1 W R 2 = f 2 AC~r 2 X R 3 = f 3 WC - r 3 Y R i = f A XB-r A Y 
R 5 = hCS - r 5 D R 6 = f 6 AD - r 6 Z R 7 = f 7 XS - r 7 Z 

(5.2) 
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Define the vector of concentration of chemical species as u = [B, C, D, S, W, X, Y, Z] 
and 



/( u > A) — [Ri, i?2, R3, R4, R5, R6, R7} , 



(5.3) 



where (-) T denotes transpose. The variations of the concentration of species on biosen- 
sor i can be expressed as [13] 

d ^=Mf(u l ,A) for t>t i: u i (t i ) = u , t = l,2,...,N. (5.4) 
where M is a 7 x 7 constant matrix obtained as fl3l. 



M = 



-1 








-1 














-1 


-1 


-1 























1 


-1 

















-1 





-1 


1 





-1 

















1 





-1 








-1 








1 


1 


























1 


1 



(5.5) 



The index i in U; = [B{ Cj D L Si Wi Xi Yi Zi] refers to biosensor i. From the reac- 
tions (5.1 ), the rate of adsorption of target molecules on the surface of biosensor i can 
be written as R(A, Uj) = f%ABi + f 2 Ad + f 6 AD l - r%Wi - r 2 X l - r§Zi. Thus the 



corresponding boundary condition in (2.6) can be expressed as 
dA 



7 



dz 



= Aq 1 u, 



T 

P , 



i = 1,2,...,JV, 



(5.6) 



z=0 



where q and p are vectors which can be defined as q = [/1 f 2 Jq 0] T and 
p = [0 r\ r 2 r 6 ] T . Applying a small alternative potential between the gold 
substrate and a reference electrode in the test solution generates a charge at the gold 
surface which causes electrons flow through ion channels 20 . The measured current 



in the external circuit is proportional to the surface average of dimer concentration. 
Denoting the average dimer concentration of biosensor i by Di(t), the observation 
equation, on biosensor i, can be written as mi(t) — Di(t) + rii(t). Physical reality 

(|5.4|) , and (|5.6|) has a non- 



demands that the derived PDE model defined by (2.2 )-( 2.5 ) 
negative solution since the solution corresponds to non-negative physical quantities. 
In the following, Theorem |5.1| states that the solution is actually non-negative. By 
proving the positivity of the solution, it can be verified that the PDE model (|2.2|)-(2.5), 



and (5.4), (5.6) is well defined. 



Theorem 5.1. Let (2.2)-(2.5), \5.A) , and (5.6) describe the dynamics ofA(t,y,z) 
and Ui(t) fori = 1, 2, . . . , N. Then, Ui(t) > for t > U and A(t, y, z) > for y G [0, 1], 
z £ [0, h], and t>0. 

Proof. The proof can be found in Appendix [Cj □ 

5.2. Illustration of the accuracy of the multi-compartment model. This 
section considers a biosensor array comprising four ICS biosensors. The aim is to 
show that the multi-compartment model of Sec |3.2| yields an excellent approximation 
to the flow dynamics. The height and width of the flow chamber are h = 0.1 mm and 
w = 2 mm. The length of each biosensor is L — 2 mm and the diffusion constant 
is equal to 7 = 10~ 6 cm 2 /s. We study the effect of varying the concentration A* at 
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'ft) 

.0.06 — 

: ..... . e (f) 



■7 



eft) 




eft) 




— eft) 
eft) 









400 600 
time [s] 
(a) 



400 600 
time [s] 
(b) 



F ig. 5 .1. Th e multi-compartment ODE model (3.5) is compared with th e PD E model 
\2.2\ -(2.5), '5.4-), and (5.6) by plotting the normalized error described in \5. 7| ) for four 
biosensors. Plot (a) shows the results for A* = 10~ 8 Mol/m? and plot (b) corresponds to 
A* — 10 -11 Mol/m 3 . The flow rate is set to 10 /iL/min. The length of the biosensors is 
L — 2 mm and their spacing is d = 1 mm. 



the inlet of the flow chamber in the range 10~ n ~ 10~ 8 Mol/m 3 and the velocity in 
the range 10 ~ 100/iL/min. The notation Mol/m 3 , throughout the paper, stands for 
mole per meter cube. The spacing between biosensors is d — 1 mm. 



The response of biosensors obtained by the multi-compartment model (3.5) is 



simulated and compared with the response of the PDE model (2.2|-(2.5|, (5.4), (5.6). 



The Comsol multi-physics simulation software is used to solve the PDE via the finite 
element method. The predefined convection and diffusion application mode in 
Comsol is used to define the governing PDE in the domain. The ODEs (5.4) on the 
boundary are defined through the weak form boundary setting. 

Since the measured output of the biosensor is a linear function of the average dimer 
concentration, define the normalized error between the ODE and PDE responses as 



(t) = A(t) - A- (*) /AW, < = i, 



(5.7) 



mo 



Ds( t) is th e av erage d ime r concentration on biosensor i, obtained by the PDE 



2.2)-(2.5), (5.4), and (5.6) and A° DE W i s the corresponding response from 
the multi-compartment model (3.5). Fig 5.1 shows the normalized error (5.7) versus 



time for two values of concentration A* = 10 11 and A* = 10 8 Mol/m 3 . It can be 
seen that the error during 1000 seconds of simulation time is less than 0.015% for 
A* = 1Q- 11 Mol/m 3 and less than 8% for A* = 10~ 8 Mol/m 3 for all the biosensors. 



Fig 5.2 shows the normalized error (5.7) for the first and second biosensor for dif- 
ferent flow rates. The concentration of target molecules is A* = 10~ n Mol/m 3 . The 
figure shows that by increasing the flow rate to 100/iL/min, the multi-compartment 



model (3.5) remains accurate within 9% error for the first and second biosensor. 



5.3. Investigating the properties of the estimator. This section compares 
the error variance of the estimate A\ (using the approximation (4.3)) with numerical 



simulations. Here, the analytical results of Sec(4]for the ICS biosensor array is com- 
pared with the corresponding simulated results. The achievable improvement in the 
estimate based on the number of biosensors is also evaluated. 



The variance cr| N of the finite-sample estimator A\ 



(4.1 ) is estimated by Monte 



Carlo simul ation s. The results are then compared with the approximate value (4.3) 
in Corollary 4.5 for verification. The standard deviation of A\ for different number 
of biosensors is shown in Table |5.1| The variance is obtained when S = 300 samples 



with sampling rate 1 sample/s are used for estimation. The actual value of the con- 
centration is A* = 10~ 8 Mol/m 3 . The signal to noise ratio, defined as the ratio of 
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F ig. 5 .2. Th e mult i- compartment ODE model (3.5) is compared wi th th e PDE model 



(2.5), , 5.4), (5.6) for different flow rates. The normalized error (5.1) is plotted for 



le first (plot (a)) and second biosensor (plot (b)). All the results correspond to A* 
10" 11 Mol/m 3 . The length of biosensors is L — 2 mm and the spacing is d — 1 mm. 



initial dimer concentration squared to the noise variance, is set to 10 dB. Table |5.1| 
shows that the estimation variance, for N = 1,2,3 biosensors, is decreasing slightly 
less than l/N as the number of biosensors increases. The reason is that di, defined 
in (4.3), is increasing for i G {1,2,3} in this case. In order to justify this result, the 



behaviour of the response g(A,t), defined in (2.8), is investigated. Fig 5.3 illustrates 



the response of ICS biosensor g(A,t) in (2.8) versus the concentration A of target 
molecules in the outer compartment. Recall that the measured output current of the 
biosensor is proportional to the average dimer concentration which decreases as the 



concentration A of target molecules increases. Fig 5.3 shows that the response curve 
is sigmoidal and saturates for low and high concentration. As a result, there exists 
a concentration A c (t) for each time instant t such that in the range A > A c (t), the 
response g(A,t) is convex in A. Therefore, the derivative dg(A,t)/dA is increasing. 
Since dg(A, t)/dA is negative, [dg(A, t)/dA] 2 is decreasing in A for A > A c (t). Define 

H(A) as H{A) = Y,k=i [9g(A,t k )/dA\ 2 . Then, there is a concentration A* such that 
in the range A > A*, H(A) is decreasing. Assume that the sampling tim e points t l,k 
for each biosensor are selected such that the time difference t l ' k — ti in (4.3) is con stant 
and equal to t k for i E {1, 2, ... , N}. Then, the sensitivity di, defined in (4.3), can 
be written as di — a 2l ~ 2 H(a l ~ 1 A*). The behaviour of H can explain the ascending 
behaviour of di in Table 



5.1 



To this end, H(A) and a H(aA) are illustrated versus 



A in Fig 5.4 From the experimental values of the parameters, a in (3.7) is obtained 
as a = 0.8173. The number of time samples in acquiring H(A) is S = 300 and the 
sampling rate is 1 sample/s. Fig 5.4 shows that there is a range [A m ,v4 n ] such that 



2 H(aA) > H(A), Ae[A m ,A n ], H{A) = ^[dg{A 1 t k )/dA\ 



(5.8) 



k=l 



Consequently, di in (4.3) is increasing as long as a 1-1 A* belongs to [A m , i"] for 
i € 1, 2, ...,N — 1. Assuming that A* e [A m , A n ], there is a positive integer N* such 
that for N < N* , di in (|4.3|) is increasing for i = 1, 2, . . . , N. Hence, using N < N* 



biosensors decreases the estimation variance to less than l/N of the variance obtained 
with a single biosensor. The expression for TV* can be obtained as 



N* = I 



log a 



■J +2, 



(5.9) 



where A m is defined in (5.8), A* is the concentration of target molecules at the inlet 



of the flow chamber, and a is defined in (3.7). In (5.9), |_ - J denotes the floor function 



which maps a real number to the largest integer which is not larger than that real 
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Table 5.1 



Comparison betwee n the simulated and approximate value ,4-3), in Corollary 4-5 for 
the variance <r| N of A\ \4-D- The simulated and analytical values of the standard deviation 
o~s,n/A* for S — 300 time samples are shown. The sampling rate is 1 sample/ s. A* — 
10" 8 Mol/m s . The signal to noise ratio, defined as the ratio of initial dimer concentration 
squared to the noise variance, is equal to 10 dB. 





<ts,n/A* 


Simulated 


Analysis (4.3) 


N=l 


0.0971 


0.0943 


N=2 


0.0611 


0.0651 


N=3 


0.044 


0.051 



A [Mol/m 3 ] 



Fig. 5.3. The response g(A,t), defined in (2.8), is insensitive to low concentrations of 
target molecules. By increasing the concentration of target molecules, the sensitivity of g(A, t) 
(the magnitude of the derivative dg(A,t)/dA) is initially increasing and then decreasing. The 
response is plotted for a time horizon of 500s. 



number. According to the value of A™ 1 in Fig 5.4 the value of N* is obtained as 
N* = 4 for A* = 10- 8 Mol/m 3 . 



It can be seen from Fig |5.3| and Fig |5.4| that increasing the number of biosensors 
beyond a certain number does not decrease the estimation variance significantly since 
the derivative dg(A, t)/dA approaches zero as target molecules in the fluid are grabbed 
by the previous biosensors in the array. 



The asymptotic behaviour of the esti mate A\ (4.1 1, using an array of N biosensors, 
for N = 1,2, 3, is illustrated in Table I 



5.2 



in (4.2), is computed for N 



where <t| N is defined in Corollary 
variance. The variance cr| N 

io- 



4.5 



The asymptotic variance of y SA\ , obtained 
1,2,3. Then the simulated value of S r (r| JV , 
is compared with the corresponding asymptotic 
is obtained by Monte Carlo simulations. The initial 
concentration is A* = 10~ 8 Mol/m 3 . In summary, it can be seen that: 

1. The multi-compartment model (3.5) predicts the response of the ICS biosen- 
sor array accurately over the range of concentrations between 10 -11 Mol/m 3 and 
10" 8 Mol/m 3 . 

2. The initial concentration A* at the inlet of a flow chamber can be estimated 



with the ICS biosensor array using the multi-compartment model. Table 5.1 shows 
that by using an array of three sensors and 300 measurement samples during the first 
300s of the response, less than 5% estimation error can be achieved. 

3. When the initial concentration A* is in a certain range and the number of 
biosensors N does not exceed a certain limit, the estimation variance is less than 1/N 
of the variance obtained by a single biosensor. This shows that there can be significant 
advantages in using an array of biosensors compared to a single biosensor. 



Appendix A. Proof of Theorem |3.1[ Divide the flow chamber into two 
segments along the z direction at z = ho where ho is selected as (3.2). Then discretize 
the y-axis at the edges of the biosensors located at y = y^i and y — yi^i, as shown 
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10 

A [MoL'm J ] 
(b) 



Fig. 5.4. The ascending behavio ur of dj in (4-3) is investigated by comparing H(A) and 



a H(aA) where H{A) is defined in (5.8). The par t of plot (a) where a H{aA) > H(A) is 



magnified in p lot (b). The interval [A"\ A 

of N* in , 5.9) is obtained as N* = 4 for the initial concentration A 



5.8) is specified for this example. The value 

8 Mol/m 3 . 



1(T 



Table_£L2 

Asymptotic behaviour of A\ defined in (4-1)'- The finite-sample estimation variance 
a N s °f Ai is obtained by Monte Carlo simulations. The concentration at the inlet of the 
flow chamber is A* = 10 _8 /m 3 . From 4-2), cr 2 r _1 is co 







^r- 1 


S = 300 


S = 1000 


S = 5000 


N = 1 


2.8273 x 10~ ib 


7.4021 x 10- lb 


2.3283 x 10" ib 


2.491 x lO" 1 ^ 


N = 2 


1.1199 x 10~ ib 


3.2628 x 10~ lb 


1.5153 x 10~ ib 


1.139 x 10~ 12 


N = 3 


5.808 x 10- iY 


1.5981 x 10~ lb 


7.2200 x 10- ib 


0.692 x 10" 12 



in Fig. 2.2 The inner compartment of biosensor i is located at y € 1/1,2) and 
z G (0, ho). The outer compartments above biosensor i — 1 and i are separated by a 
middle compartment which is located in the range y € (j/i-1,2, Vi,x) an d z € (ho,h). 



Non-dimcnsionalizing (2.2) yields 

, d 2 a d 2 a eZ(l - Z) da 



, da 



8t dZ 2 dY 2 p dY' 1 

a = A/ A*, r = j/h 2 t, Z = z/h, Y = y/L 

The Peclet number \21 in this problem is 

P e (Z)=eZ(l-Z)p-\ 



L/h, p 



7 
4hv 



(A.l) 



(A.2) 



which is used as a criterion to identify the dominance of advection or diffusion in the 
transport of particles in the flow. Obtain the PDE (|A.1[) at Z = e /h as 



da 



Z=€ /h 



d 2 a 
dZ 2 



Z=t a /h 



a a 
W 2 



1 



£11 

h , 



Z=e /h 



da 
dY 



(A.3) 



Z=e /h 



Assume that eo = 0(j ). Therefore, the Peclet number (A.2) at Z = eo/h is P e (^) = 
0(7) and the advection term on the right hand side of (A.3 1 is negligible comparing 
to the diffusion term. Thus, (A.3) at Z — e /h is rewritten as 

, da 



dr 



Z=e /h 



d 2 a 
dZ 2 



Z=€ /h 



era 

W 2 



0(7), 



(A.4) 



Z=e /h 
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Hence, the dimensional PDE (2.2) at z = eo can be expressed as 

0(7 2 ), 



OA 




d 2 A 




d 2 A 




~dt 


z=e 


^ dz 2 


z=t 


^ dy 2 


z=e 



(A.5) 



From the selection of ho in (3.2 ), the magnitude of ho is of order 0(7 1 / 3 ). Assume 
that the value of eo is selected such that < eo < ho- The second order derivative in 
z direction in (A.5) can be expressed as 



d 2 A 


1 


' dA 




dA 




Ih 2 


*=e ~ h ° 


dz 


z=h 


dz 


z=Q 



+ O(h ) 



(A.6) 



The derivative in (A.6) can be expressed as 



OA 

dz 



z=h 



— [A(t, y, h + e ) - A(t, y, e )] + O(h ). 
h 



(A.7) 



Using the boundary condition (2.6), the derivative §7| z _ m (A.6) on sensor i is 



dA 

dz 



1 



~R(A(t,y,e ) + O(e ),Ui) , y € (y itl ,yi, 2 )- 

z=Q 7 



Substituting (|A_7j) and ( |A.8| in ( |A.6| yields 
d 2 A 



dz 2 



1 

h 
1 

hoi 



1 

h 



[A(t, y, h + e ) - A(t, y, e )] + O(h ) 



(A, 



(A.9) 



R (A(t, y, e ) + O(e ), u*) + O(h ). y £ {y iA , j/ i>2 ) 



From (A.9), the PDE (A.5) for y £ (j/i, 1,2/4,2) can be expressed as 
dA 7 



h 



dt 



h 



[A(t, y, ho + e ) - A(t, y, e )] - R (A(t, y, e ), u,) + O(e ) (A.10) 



+ lO{ho) + Ki fl . 



d 2 A 

dy 



+ iO(h 2 ) + /i O( 7 2 ), V € {Vi,i,Vi,2) 



The exponent of the concentration A in the polynomial R(A, u) can be equal or 
greater than one. Thus, the maximum order of approximation error in R(A, u) for 
an approximate A occurs when the exponent of A is one. This case is considered in 
( |A.10[ ) and R(A(t,y,e ) + O(e ),Uj) is replaced with R (A(t, y, e ), u») + O(e ). The 
value of e is selected as e = 0( 7 2 ). From (E§, h = 0(7 1/3 ). Thus, we have 
JL = 0( 7 2 / 3 ), -fh = 0( 7 4/3 ), jhl = 0( 7 5 / 3 ), and h Q ~/ 2 = 0( 7 7 / 3 ). Thus, the sum of 

the last five terms on the r ight h and side of (A. 10) is of order 0(7 4 / 3 ). The firs t term 
on the right hand side of ( A. 10 1 is o f order Opp' 3 ) and the left hand side of (A. 10) 
is of order 0(7 1 / 3 ). Equation (A. 10) for y € (2/^1,^,2) can be rewritten as 



ho 



dA 
~dt 



[A(t, y, h + e ) - A(t, y, e )] - R {A(t, y, e ), u 4 ) + 0( 7 4/3 ), (A.ll) 



where ho is defined in (3.2). 
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The rest of the proof focuses on obtaining an expression for A(t, y, ho + eo) in 
(A. Ill for y £ {Vix,Vi,2)- Consider (A.ll for Z £ (h /h,l — h /h). It can be seen 



from the selection of ho in (3.2) that the Peclet number (A. 2) has a large value for 



a small 7 and can be written as P e (Z ) — 1/A where A at its maximum is of order 
0( 7 2 / 3 ) at Z = ho/h. The PDE |Al} for Z £ (h /h, 1 - h /h) can be expressed as 



, da 
dr 



Ae 2 ^ = Ae 2 



da 2 da 2 
dZ 2 + W 2 



da 
<9Y' 



A = P~\Z), Z e (ho/h, 1 - ho/h). (A.12) 



Consider (A.12) for t £ (t»-i,tt) where U is defined in (3.4). For t £ the 



part of the flow chamber, in the range y £ (2/1-1,2,2/1,2), is partially occupied with 
the fluid. Therefore, da/dY in ( A.12 ) has a significant value such that A -1 is 



of order 0(X~ 1 ) where A, in (A.12), has a small value. Thus, comparing (A. 4) and 
( A.12[ ) concludes that for t £ (£t-i>*t) the variations of y, z), for z £ {ho,h — 
ho) and y G (2/4-1,2,2/4,2), is fast whereas A(t, y, eo) is slowly varying for y £ (0,1). 
Therefore, assumption (2) can be justified. By ignoring the diffusion along y and z 
direction versus the advection term in (A.12), the dimensionalized form of ( A.12[ ) for 
t £ (tj_i, t{), y £ (2/1-1,2, 2/2,2), and z £ (ho, h — ho) can be written as 

dA dA 

— = -v(z)— + 0(7), z £ (h , h - h ), y £ [Vi-x, 2,2/2,2), t £ (t i ^ 1 ,t i ). (A. 13) 
The above PDE solution has the following initial condition. 



A(ti-t,y,z) = 0, 2/ G (2/4-1,2,2/2,2), z £ (h , h - h Q ). 



(A.14) 



In order to find the solution of (A.13), a boundary condition is required. By applying 
the divergence theorem in the outer compartment above biosensor i — 1, a boundary 
condition is obtained at 2/i-i.2- 

Determining a boundary condition for (A.13) at y = 2/; i.2 : Consider the 



general form of the advection-diffusion equation in the flow chamber (2.1 ) as 



dA 
~dt 



(7VA 



Av 



ye (0,1), z£(0,h) 



(A.15) 



where v is the velocity vector of the fluid, V • (•) represents the di vergen ce operator, 
and V denotes gradient. Applying the divergence theorem [22] to (A.15) results in 



dA 
~dt 



dy dz 



(jVA 



Aw 



dl(y,z) 



(A.16) 



where the left hand side is a surface integral over the bounded domain ft. The right 



side of (A.16) is a line integral over the boundary of ft. The vector n is the outward 
pointing unit normal field of the boundary dft. Here, the domain f2 is selected to be 
the rectangular region defined as ft = {(y,z)\z £ (ho,h),y £ (2/1-1,1, 2/4-1,2)}- The 
domain f2 comprises the outer compartment of biosensor i — 1. Assuming that the 



multi-compartment model (3.5) holds for biosensor i— 1, the maximum time derivative 
of the concentration in the outer compartment above biosensor z — 1 is 0( 7 ). Thus, 



the left hand side of (A.16) is of order 0(7). The line integral, on the right hand side, 
over the upper edge of Q at z = h is zero due to insulation boundary condition in 
(2.5). The middle compartment between sensors i — 1 and z — 2 is also in equilibrium 
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with concentration Ai-x + 0(7). Hence, the line integral over the left edge of fi in 
(A. 16 1, at y — is obtained as 



(jVA - A-v^j ■ {-ydz) = / v{z) dz + 0(7). 



(A.17) 



Denoting the concentration in the inner compartment of biosensor i — 1 by ai„i(i, y), 



the derivative j^\ z - hg can be expressed as — 1 ^ l( - t ' t '- ) + O(ho). Using this expres- 
sion, the part of the line integral in (A.16), over the bottom edge of fl at z = ho, can 
be written as 



(jVA - Av) ■ (-My) = ^ (o<_i(t) - M-i) + 0( 7 4/3 )- 



h 



(A.18) 



Here, aj_i(f) is the average of <2i_i(i, y) over the length of the biosensor. By sub- 
stituting ( A.17| ) and (A.18 1 in (A.16), the line integral over the right edge of Q at 
y = 2/1-1,2 is evaluated as 



h / OA 

J-Qy ~ Av(z) ) dz = -f(t), y = u, 1 



(A.19) 



where h 

f(t)=A^ 1 (h-h )v 1 + ^ r (a^ 1 (t)-A i _ 1 ) + O( 1 ), v 1 = J ^ i l. (A.20) 

ho h — nv 

Since the concentration everywhere in Q is equal to A^_i + 0(7), the derivative in z 
direction along the left edge of is of order 0(7). Therefore, (A.19) results in 

OA 



7 



dy 



Avi = - 



/(*) 



, V = Vi-1,2, z G (h ,h). 



(A.21) 



h — ho 

where f(t) and v\ are defined in (A.20). Note that the term 7^ cannot be considered 
to be of order 0(7) since ^ at t = U_i and y — j/i-1,2 can be of order 0(7 _1 ) or 
larger. Equation (A.21 ) can be considered as a boundary condition for the transport 
equation (A.13). 



The solution of the transport equation (A.13) with the initial condition (A. 14) 
and the boundary condition (A.21) for t e [U-i, ti) can be found as [23] 

A(*,y,z) = h(t - U-i,y - Vi-i,2,z) + 0(7)) v 6 [y»-i,2,s/»,2), z e (h a ,h- h ), (A. 22) 

where, 



h(t,y,z) 




exp 



,t-- 



j,, exp(^r)/(r + V 1 )dr 



^(y_l,(z)i)j/ ^ 

y > > 



v f(y-v(z)t) 



Using (A.20) and integration by parts obtains h(t,y,z) in (A. 22), for y < as 



h(t,y,z) = aAi-i + 



f(U- 



jL 



h (h - h )v l 



exp 



'T 



(h-h )vi vi(h-h )Jo 



exp 



u(z) 

<f/( r + tj-i) 
dr 



(y ~ v(z)t) 



dr 



(A.23) 

0(7), 
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where a is defined in p.7[). Consider (A.23) for y < tv(z) — ftj 1 / 3 where /3 is a pos 



itive constant. In ( |A.23| , d ^j^ can be replaced with ^ da ' d ^ T>> according to (|A.20|) 
According to assumption (1). 



daj_ i (t) 
dt 



dai- 1 {t) 
dt 



> for ti-i < t < ti. Thus, the maximum 



is obtained as ^Ai-i from (3.5) since R(A,u) > 0. Considering 



derivative 

assumption (3) and y < tv(z) — /37 1 / 3 , the maximum magnitude of the third term 
on the right hand side of (A.23) is of order 0(7). From the above explanation, the 
maximum change in <2j_i, for t € [t^_i,tj), is Av4.i_i(ii — fi-i). Therefore, according 



to assumption (3), the second term on the right hand side of (A.23) is also of order 
0(7). Consequently, it can be expressed that 

h(t, y, z) = aAi-i + 0(7), V < Hz) - /? 7 1/3 - 



Considering (A.22) and the above equation, the concentration A(t,y,z) inside the 



outer compartment of biosensor i and the previous middle compartment at t = i, 
(when the flow reaches the far end of biosensor i) can be obtained as 



A(ti,y,z) = aA t -i + 0(7), y G [yi-1,2, Vifi), z € (h ,h - h ), 



(A.24) 



where a is defined in (3.7). It is shown at the end of this section that for t > ti the 



variations of derivative"^ for y e (yi-1,2, 2/1,2) an d z £ (/io, ft — ho) is of order 0(7). 
Therefore, there exists a time t* such that for t € (ti, t*), ( |A.4[ ) holds. The complete 
proof of (A. 4) can be found at the end of this section. The value of A(t, y, ho + eo) in 
(A.ll), for t £ (ti,t*), can be obtained as A(t, y, ho + £0) = Ai + 0(7) according to 



(A. 4) and (A.24). Hence, A(t,y,e ) in (A.ll) forte (U,t*) is governed by 



h 



dA 

~dt 



= ~[A- A(t, y, eo)} - R {A(t, y, e ), u t (t, y)) + 0( 7 4/3 ), (A.25) 
y e {yi,i,yi, 2 ), tefat*), A(t u y,e ) = 0. 



where Ai is obtained in (3.7). Considering assumption (2) and eo = CK7 2 ), (2.7) can 
be expressed as 



du t (t,y) 
dt 



G(ui(t,y),A(t,y,e ) + 0(^ 2 )), te(U,t*), y e [y t ,i, y il2 ] , u, (/,, y) = u . 



For an arbitrary value yi for y, (A.25 ) together with the above equation will represent 



an ODE system whose solution includes expressions in terms of t and independent of 
yi for A(t,yi,eo) and Ui(t,yi). In fact, A(t,y, eo) and Hi(t,y) at every point y has 
the same dynamics in time. In order to derive (3.5), A(t,y,e ) and Uj(t, y) in ( A.25[ ) 
and the above equation are replaced with their average values in the range (yi,i, yi,2)- 
The spatial average of A(t,y,eo) in (yi,i,y-i,2) is denoted by a,(f). 

To complete the proof, it needs to be shown that the multi-compartment model 



(3.5) holds for the first biosensor which is straightforward. The concentration in the 



outer compartment of the first biosensor achieves equilibrium fast. For t > ti, it is 
equal to the concentration A* at the inlet of the flow chamber. 



Proof of ( A.4): From assumption (1), A(t, y, z) is an increasing concave function 
for t > U and y £ (yi-1,2, Vi,2 
y € (yi-1,2, Vi, 2 )- By showing that 



in t for t > ti and y G (yi-1,2, yi, 2)- Thus ^ > is decreasing in t for t > U and 



&4(*i,y,z) 
dt 



0(7), y e (yi-1,2, y», 2 - 7 1/3 ], z e (/i , ft - fto), 



(A.26) 



BIOSENSOR ARRAYS FOR ESTIMATING MOLECULAR CONCENTRATION IN FLUID FLO\f 9 



it can be concluded that ^ 



= 0(7) for all t > ti and thus (A.4) is obtained. For 
the derivation of (A. 26), the value of -j^A (t, y, z) for y g (2/2-1.2,2/1.2) at t = ti can be 
obtained from QA.22p , ( |AT23| ), and ( |A.20| as 

dA{ti,y,z) _ jL ddi-! u Ay 



v{z)- 



exp 



dy 

—(Ay ■ 
7 

1 L 



hoVi(h — ho) dt 
v{z)v\k 



v(z)At) 



h vi(h - 
where Ay = 



h ) 



exp 



—(Ay 
7 



v(z)At) 



v(z) 

v(z)f(U-i) 
j(h - ho) 

At 4*. 



(A.27) 



7L doi-i^j-x) 



exp 



/i fi(/i - h ) 
v\v(z)t d 2 di-\ 



7 



dt 2 



dt 



V - J/i-1,2 and A/ = 



The first term on the right hand side 



of (A.27) is of order 0(7). The second term converges to zero faster that 0(7) for 
V < Vi.2 ~ /?7 1 ^ 3 - Here, /3 is a positive constant. Assuming that <Zj_i(t) is con- 

cave an d 1 < 0, the maximum magnitude of the integral on the right hand 
side of ( A.27| is obtained by setting the exponential term in the integral equal to 



vivjzj 
7 



(At — Ay/v(z)). Therefore, the maximum magnitude of the last term on 



exp 

the right hand side of (A.27) is of order 0(7). Thus, 



v{z) 



dA(t u y,z) 
dy 



0(7), y G (2/»-l,2iS/i,2 -7 1/3 ]j z £ (h ,h- h ) 



By substituting the above equation into (2.2), (A.26) is obtained. 

Appendix B. The proof of Theorem |4.4| on the asymptotic properties 
of the estimator. For the proof of part (1), it is required to show that gi(Ai,t) 
in (4.1) is continuous in A±. In the multi-compartment model (3.5), and di are 
continuous in A4 since R and G are continuous functions in u^, di, and Ai and Lip- 



schitz in and di with Lipschitz constant independent of Ai 24 . From (2 
gi(A\,t) = F(\ii). Since F models an electrical response as a measurable physical 
quantity, it is continuous in U;. Thus, gi(A\,t) is continuous in u^. In our case 
study in SecjEJ F is equal to one of the elements of u^. On the other hand, from 
the recursions in (3.7), Aj is continuous in A±. The continuity of gi(A\,t) in u^, 
Ui in Ai, and Ai in A\ concludes that gi(A\,t) is also continuous in A\. Relating 
the current estimation problem to the one stated in Theorem |4.2[ f t in the theorem 



,g 2 (A 1 ,t l > k ),...,g N (A 1 ,t l > k )] . The noise vec- 

S are independent identically distributed with 
it can be shown that the 



4.2 



corresponds to the vector ( Ai , t 
tors \n\ , nf , • • • , n%] T for k = 

covariance a 2 I. Dealing with condition (b) in Theorem 
tail cross product of g%(A\, t*' ) exists. It is required to show that for each A\,A 2 > 0, 
S^ 1 X)fe=i ffi(^i) t l ' k )9i{^2, t l ' k ) converges uniformly to a specific function of Ai and 
A2 when S, t l,s — > 00. When t l ' S — > 00, the concentration of chemical species on each 
biosensor eventually achieves equilibrium. So does the response gi which is a function 
of the concentration of chemical species. Thus, gi(A,t 1 ' ) — > g e ^(A) where g e ^(A) 
is the steady state value of the response. Assuming that the limit function g e ^(A) 
is continuous, it can be concluded that gi(A,t l ' S ) — > g e ^(A) uniformly as S — > 00. 
Therefore, S' 1 £fe=i 9l {A u t^ k ) gi (A 2 , t> k ) converges uniformly to g e ,i(Ai)g e ^(A 2 ). 
In order to satisfy assumption (c) , the existence of corresponding tail cross products 
can be similarly proved. 

There is a one-to-one correspondence between the value of concentration Ai and 
the concentration of chemical species at equilibrium. Since there is also a one-to-one 
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Fig. C.l. Discretizing the flow chamber along the length and height of chamber. The 
length and width are respectively discretized to P + 1 and Q + 1 steps 



correspondence between the response and the concentration of chemical species on 
the biosensor, gi(Ai, t l,k ) is a bijection function of A\. Therefore, the only value of 
A 1 which minimizes Q{A X ) = lim s ^oc S' 1 £f =1 Ylti (g t (A* , t^ k ) - 9i (Ai, ?'•*)) 2 
Ai = A* . Thus, condition (b) in Theorem ■ 



In our case, a(9o) in Theorem 4.3 corresponds to 



4.2 



giiAi,?'-)) is 
is satisfied and A\ is strongly consistent. 



r = lims-^ S- 1 Eti Ef=i {d 9i {A*,t' k )/dA 



(B.l) 



where dgi(A* , t l > k )/dA 1 is the value of dgt(Ai, t l ' k )/dA 1 at Ax — A*. We assume that 
all biosensors have identical parameters and binding properties. Thus, according to 



(3.5), there is a common functional relation between the response of each biosensor 



and the concentration in its outer compartment expressed as 

g i (A 1 ,t)=g(A i ,t-t i ) for i = l,...,N 



(B.2) 



Here, g(Ai,t) denotes the response of any biosensor when the concentration in its 
outer compartment is A4. The response g[A,t) commences at t = 0. Recall the time 
shift ti in ( |B.2[), i s the response delay of biosensor i. According to the recursion (3.7 1 
in Theorem |3.1| (B.2) can be expressed as gi(Ai,t) — g{a % ~ 1 A%,t — ti). Therefore, 
r in (B.l) can be expressed by (4.2). Regarding Theorem 4.3 the estimation error 



A\ — A* is asymptotically normal. 

Appendix C. Proof of Theorem |5.1[ In order to prove the positivity, we 
rely on the fact that the solution of the spatially discretized problem converges to the 
solution of the original PDE problem when the spatial discretization step is sufficiently 
small. Thus we can conclude the positivity of the solution of PDE by proving that 
the solution of the discretized system is positive. We use the following theorem to 
prove the positivity of the discretized system: 

Theorem C.l. Suppose that F(t,v) is continuous and satisfies a Lipschitz con- 
dition with respect to v. Then the system w'(t) — F(t,w(t)) is positive iff for any 
vector v € K m and all i — 1, to and t > 0, 



v > 0, v, = => Fi(t,v) > 



(C.l) 



By positivity of the system w'(t) — F(t,w(t)) we mean that for any initial value 
iu(0) > the solution is w(t) > for all t > 0. The requirement of a Lipschitz 
condition can be slightly relaxed. It is sufficient that the initial value problem has a 
unique solution for any w(0) > |25|. The proof of this theorem can be found in 21 



Consider the General PDE model described in Sec. [2] It can be proved that the 
solution of this model exists and is unique 23 . We discretize the flow chamber in y 
and z directions to P + 1 and Q + 1 steps respectively as shown in Fig. |C.1| The inlet 
locates at y = which is equivalent in the discretized system to i — and the outlet is 
at i — P. The array of sensors is located at j = and the upper face of the chamber 
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is at j = Q. By applying first and second order approximations, the discretized form 
of the advection-diffusion equation inside the chamber can be written as 



, .n Ai+ij — Aij Ai j-i — 2Aij + Ai A+i 
v(j) ^ -+7- 



Ay 



Az 2 



(C.2) 



for i = l,...,P-l,j = l,...,Q-l 



Where the discretization steps in y and z direction are respectively referred by Ay > 
and Az > 0. We need to consider Aij at the boundaries separately through the 
boundary conditions. The boundary condition at y — simply converts to 



.4 



Ai, j = 0,l,...,Q. 



The discretized form of boundary condition at y = I can be written as 
Ap,j - A P -i t j 



Ay 



= 0, j =0,1,. 



similarly we have the insulation condition on the upper side of chamber; 

- l; "- l - = o i = 0,l,...,P. 



Az 



(C.3) 



(C.4) 



(C.5) 



At z — 0, at points which are not located on sensors, the boundary condition can be 
written as 



A M - A lfi 
Az 



= 0forie {0,1,...,P}-U^ =1 T4 



(C.6) 



where Vk is defined as the set of indices of the discretized points which are located on 
the fc-th sensor; 



V k = {i G {0, 1, P} I zAy G [y kl i, y k>2 ]} 



The boundary condition on sensors, according to (5.6), can be written as 
A41 - A ifi 



7" 



A,: 



= A ifl q f ul - p f ul, i G V k , for i = 1, N. 



(C.7) 



(C.8) 



where \i l k = u k (t,iAy) for i G Vk- The system of ODEs in (5.4) remains unchanged 
for the discretized domain. At each point, the time derivative of the chemical concen- 
tration vector on fc-th sensor can be written as 



ul = M/(i4, A,o) *£V kl for i = 1, N. 



(C.9) 



For proving the positivity (short for 'non-negativity preserving') of the solution, we 
need to show that Ai t j > for i = 0, P and j = 0, Q and u| > for i G V k where 
fc = 1,2,..., A. Using Theorem C.l we select v to be a vector containing Aij for 
i = 1, P — 1 and j = 1, Q — 1 and u^. for i £ V k and k = 1, A as its elements. 
v can be represented as 



3 Jie{i,...,f-i},je{i,...,Q-i}: l u k Uev h ,ke{i,...,N} 



(CIO) 



Apparently, F(t,v) contains their derivatives given by (C.2) and (C.9). Consid- 
ering Theorem C.l, it is required to show that the statement in (C.l) holds for this 
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case. To this end, we split the elements of v into several sets and investigate (C.ll 
for each one separately. First set is denoted by v\ and contains the elements Aj j 
for i = 2, ...,P - 2 and j = 2, Q - 2 such that vi = {A itj \i G {2, ...,P - 2},j G 
{2, Q — 2}}. We should show for Ajj € t>i that AJ ■ is not negative when Ajj = 



and v > 0. Considering these assumptions and according to (C.2), A\ j is obtained as 



i = 2,...,P 



~-v(j)- 



A 



A 



Ay 2 
2,i = 2,. 



Ay 



(C.ll) 



.,Q-2 



and 



We have i € {2, ...,P-2} and j G {2, Q-2}. Therefore, Ai-ij, A i+1J , A, 
Aj ; j + i belong to the vector v and are positive according to the assumption in (C.l). 



It is easy to show that for Ay < j/v, A[ j in (C.ll ) is non-negative. The second part 
of v is denoted by V2 and comprises the values oi A± t j for j = 2, Q — 2; 



v 2 = {A ld \j€{2,...,Q-2}} 



(C.12) 



According to (C.2 1, A[ ■ is written in terms of Aqj, A2J, Aij_i, and Aij+i. For 
the values of j e {2, Q — 2}, all these elements, except Aoj, belong to the vector u 
which is assumed to be non-negative. A 0j - does not belong to v but it is non-negative 



due to the boundary condition described in (C.3). Similarly, it can thus be shown 



that for Ay < j/v, the time derivative of the second part of v which is given in ( C.12 1 
is non-negative. The third part of v contains Ap—ij for j = 2, Q — 2 and is denoted 
by t>3 = {Ap-\ t j\j € {2,...,Q — 2}}. In a similar way, it can be shown that elements of 
V3 are non-negative. To prove this, we need to show that Apj which does not belong 
to v is non-negative. According to the boundary condition in (C.4|, it is obvious that 



A 



Ap-i. 



(C.13) 



Hence, Apj = Aj 
as 



On the other hand, the assumption in (C.l ) implies that Ap_ij should be set to zero. 

■j — Ap_i t j — according to (C.13). Consequently, A' P _ 1 _j can be written 



A' 



A P _ 



p-1,3 = 7- 



P-2,j 



+ 7 



^p-i,i-i + Ap-x,j+\ 

A^ 



(C.14) 



Since all the terms on the right-hand side of (C.14) is non-negative, A' P1 ■ for j = 
2, Q — 2 is also non-negative. The fourth part ot v is considered as V4 = {A^iji G 
{2, P — 2}}. Its time derivative can be written as 



A; 



A 



A U - 7 Ay 2 
for i = 2,...,P-2 



i+1,1 



«(1) 



A 



Ay 



7- 



Aj,o + Aj, 2 
Az 2 



(C.15) 



In (C.15), Aj,o is the only element that is not included in v and its positivity should 



be investigated. According to the boundary conditions in (C.6) and (C.8), A Jj0 can 
be written as 



Ai 



N 



A ~0£t ^ fOT k = l 

Ai,i iG{0,l,-,P}-Uf =1 7 fc 



(C.16) 



where 14 is defined in (C. 7 (According to our assumption, Aj i = 0. Besides, consid- 
ering the assumption that w > 0, we have u| > 0. Therefore, (C.16) shows that A^o 
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is not negative. It can be shown that for the same small value of Ay < 7/11, A[ x > 
for« = 2,...,P-2. 

The fifth part of v to investigate is described as U5 = {Ai : Q-i\i £ {2, P — 2}}. 
According to (C.ll and (C.2|, We obtain 



A' 



Ai,Q-2 + Ai^Q 



Ai-\,Q-1 + A 



i+l,Q-l 



v(Q-iy- 



A 



i+l,Q-l 



Ay 



(C.17) 



+ 7 



Az 2 



for i = 2,...,P 



In (C.17 1, the only term that does not belong to the vector v is A^q. According to 
the insulation boundary condition in (C.5), Ai Q is equal to A it Q_i which is assumed 
to be zero. Therefore, (C.17) is reduced to 



A' 



7 



A, 



+i,Q-i 



Ay 2 

7 ^q^ for i = 2 p_ 2 
Az z 



v(Q-l) 



A, 



Ay 



(C.18) 



and according to the previous explanations, it is easy to realize from (C.18) that 
A\ q_ 1 is not negative for Ay < 7/v. The remaining elements of v for checking the 
statement in (C.l) about, are Ai t ±, Ax,q_i, Ap_li, Ap_i,q_i, and uji. for i £ Vu and 
k = l,...,N. 

For A\ 1, it can be written that 



., A 0l i + A 2) i , 1N ^2,i , ^1,0 + ^1,2 



Ay 2 



Ay 



Az 2 



(C.19) 



^2,1 and Ai } 2 belong to t; and are therefore non-negative. Regarding (C.3), Ao,i 
> 0. A\ is achieved by either (C.6) or (C.8) as 



^1,0 



A14 



1 e V k for jfc = 1, 

1 i ^ =1 v k 



(C.20) 



It is concluded from (C.20) that Ai t o is not negative. Consequently according to 
JCTTgb , > for Ay < 7/zL 

Equation (C.2) for q_i is written as 

- ^ Q - i A + y ^- Q - 1 -v(Q-D^+ 7 Ai ^:^ (C.21) 



All terms on the right side of ( C.21 ) other than Ao.q-i and Ai.q are included in v and 
are non-negative. According to the boundary condition at inlet in (C.3), Aq^q—i = 
A\ > 0. Aiq is equal to Axq-\ according to (C.5) which is assumed to be zero. 
Consequently A\ q_ 1 can be proved to be non-negative for the same small value of 
Ay. 

Considering (C.2) for i = P — 1 and j = 1 and setting Ap-\_\ to zero, we can 
state that 



Ap-i tl = 7 
+ 7 



Ap-2,1 + Api 



Ay 2 

Ap-1,0 + Ap-1,2 
Az 2 



v(l) 



A PA 
Ay 



(C.22) 
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Ap t \ can be shown to be non-negative according to (C.4). According to (C.16) and 
the assumption that Ap_ u = 0, Ap_x,o is greater than or equal to zero. The other 
elements are all included in v and are positive according to assumption in (C.l). 



Therefore, it can be easily seen from ( C.22 1 that A' P _ 1 1 
of Ay. 

For ip-i.Q-i, we rewrite (C.2) for i = P — 1 and j = Q — 1 as 



> for the same small value 



A 



p-i, 



7- 



A P _ 2 n-l + A 



P,Q-1 



Ay 2 



v(Q 



1) t 1 h 7 



.4 



p-i,. 



Ay 



According to boundary conditions in (C.4) and (C.5), Ap.Q_i = 
Ap-i,Q = Ap_x t Q_i = 0. Therefore, (C.23I can be rewritten as 



Az 2 



Ap-l,Q-\ 



(C.23) 



and 



A' 



A P _ 



P-l.Q-l 



7" 



P-2.Q-1 

Ay 2 



+ 7 



Ap-1,Q-1 

Az 2 



which shows that A' P _ 1 q_ 1 is positive since Ap-2,Q~i and j4p_i_q_2 belong to v and 
are positive. Reminding the definition of v in (C.10), we should check the validity 
of (C.l) for the remaining elements u\ when i £ 14 for k = 1,...,N. u k is the 
concentration vector at point i on fc-th sensor. Since the proof of (C.l) for each 



element of this vector is the same and independent of the others, we refer to each ui 
by u for simplicity, u has eight elements and can be presented as u = [m u 2 ...Ms]. To 



complete the validity of (C.l I, we should prove that when v > and Uj = 0, the time 
derivative u/ is not negative for j = 1, ...,8. According to (C.9), u/ can be written 



as 



Uj 



MJ(u,A lfi ) 



(C.24) 



where Mj is the j-th row of M. According to (C.16 1, Ai.o on the sensor can be written 
as 



Ain 



We have already proved that A^\ is positive for i = 1, P — 1. It is easy to show 
that Aq i and Ap i are positive. Consequently, the assumption that u > results 
in Ai t o > 0. On the other hand, by considering the elements of M in (5.5) and the 



definition of /(u, A) given in (|5.3| and (5.2), it can be easily shown that u'- > when 
we set u > and Uj — 0. From chemical point of view, the only negative terms 
in (C.24) which reduce the amount of u, correspond to the ones that have ztj as a 



tiplicative factor. This is due to the fact that the backward reaction only happens 
when Uj is not zero. However, in our assumptions Uj is zero and the negative terms 



are therefore omitted from the right side of (C.24) 



The positivity of the variables Ajj for i = and j = 1, Q— 1 can be easily shown 
since Ao t j is equal to A\ at the inlet and therefore Aqj > for j = 0, Q. According 
to the boundary condition in (C.4), Apj = Ap—ij and since we proved that Ap—i j 
for j = 1, Q — 1 is positive, we can conclude that Apj > for j = 1, Q — 1. 
For j — and i — 1,...,P — 1 we can prove the positivity by referring to (C.16) 
because we have proved that Ai t i > for i — 1,...,P — 1 and u l k > for i E Vk 
and k = 1, N. For j = and i = P, we can use the boundary condition in (C.4) 
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and write Ap_o = -Ap-1.0 an d since we just proved that Ap_x,o > 0, it can be easily 
seen that Ap o > 0. For j = Q and i = 1, P we can use the insulation boundary 
condition given in (C.5) and write A+q — A^q_i. We have already proved that 
Ai y Q-i > for i = 1, .... P — 1 using Theorem C.l The proof of positivity of Ap_Q_i 
is explained earlier in this paragraph. The proof of positivity of the solution of the 
PDE model of Sec. [2] is now complete. 
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